% Compute the global solution using a paramerized expectations algorithm
% Dynamic bank run model
% Bank leverage choice only
% Capital depreciation shock

close all
warning off

load('rbc')

%% Generate random shocks

options = optimset;
T = 10000;               % simulation periods
t_drop = 1000;

rng(1)                   % same random shock
shocks = normrnd(0,1,T+1,1);  % random shock
shocks = zeros(1,T+1);

kv   = zeros(T+1,1);
rbv  = zeros(T+1,1);
nv   = zeros(T+1,1);
xiv  = zeros(T+1,1);
hv   = zeros(T+1,1);
cv   = zeros(T+1,1);
yv   = zeros(T+1,1);
iv   = zeros(T+1,1);
lv   = zeros(T+1,1);
pibv = zeros(T+1,1);
runv = zeros(T+1,1);
rhseev = zeros(T+1,1);
rkhat  = zeros(T+1,1);
rkhats = zeros(T+1,1);
prv  = zeros(T+1,1);
av   = zeros(T+1,1);
erk  = zeros(T+1,1);
x    = zeros(T+1,npol);
x_nonlog    = zeros(T+1,npol);
x_allnonlog = zeros(T+1,npol);

%% Simulation

kv(1)  = k_ss;
av(1)  = a_ss;

    for t=2:T+1
        av(t) = rho_a*av(t-1) + sigma_a*shocks(t);
        hv(t)  = ((1-alpha)/psi*exp(av(t))*kv(t-1)^(alpha))^(1/(alpha+1/nu));
        yv(t)  = exp(av(t)) * kv(t-1)^(alpha)*hv(t)^(1-alpha);

        x0 = [log(kv(t-1)) av(t)].^expmx;
        x(t,:)  = prod(x0,2)';
        
        rhseev(t) = exp(x(t,:)*solved_eta1_rhsee);
        cv(t) = rhseev(t)^(-1) + psi*hv(t)^(1+1/nu)/(1+1/nu);
        iv(t) = yv(t) - cv(t);
        kv(t) = (1-delta)*kv(t-1) + iv(t);

        rbv(t) = 1 - delta + alpha*yv(t)/kv(t-1);

    end
    
stss_y_rbc  = yv(5000);
stss_c_rbc  = cv(5000);
stss_h_rbc  = hv(5000);
stss_i_rbc  = iv(5000);
stss_k_rbc  = kv(5000);
stss_rb_rbc = rbv(5000);

return
